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The formalism for exactly calculating the retarded and advanced Green's functions of strongly 
correlated lattice models in a uniform electric field is derived within dynamical mean-field theory. 
To illustrate the method, we solve for the nonequilibrium density of states of the Hubbard model in 
both the metallic and Mott insulating phases at half-filling (with an arbitrary strength electric field) 
by employing the numerical renormalization group as the impurity solver. This general approach 
can be applied to any strongly correlated lattice model in the limit of large dimensions. 
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Introduction. The many-body formalism for nonequi- 
librium problems was formulated independently by 
KadanofF and Baym [ij and Keldysh Q in the 1960s. 
One of the main applications of that work was to deter- 
mine the nonlinear transport properties of strongly cor- 
related materials. Recently there has been a significant 
emphasis placed on examining small open systems (quan- 
tum dots attached to leads) within the Meir-Wingreen Q 
generalization of the Kadanoff-Baym-Keldysh approach 
and there has been progress in applying Bethe ansatz 
techniques to this problem [4], numerical renormaliza- 
tion group techniques within a scattering state formal- 
ism j5l| , and Hirsch-Fye quantum Monte Carlo techniques 
by mapping to an effective imaginary time formalism Ql . 
That work is timely because the size of local electric fields 
placed over nanoscale electronics is enormous and non- 
linear quantum effects are likely to be critical for under- 
standing how those systems behave. In this work, our 
focus is on larger systems (bulk materials) placed under 
large electric fields, which serves as a counterpart (top- 
down) approach to the problem, and could have direct 
application in ultracold atomic systems placed in optical 
lattices and driven into nonequilibrium by accelerating 
the lattice through space (by detuning the frequencies of 
the counterpropagating lasers, or simply from the force 
of gravity). 

There has also been much effort applied to under- 
standing the original Kadanoff-Baym-Keldysh formalism 
and how to approximately solve the resulting equations. 
The generalized Kadanoff-Baym approximation j7j and 
the reconstruction theorem for the lesser Green's func- 
tion Q have provided much insight into the way quan- 
tum systems relax and ultimately reach a steady state. 
But there remains no exact solutions for strongly corre- 



lated bulk systems placed in large electric fields in the 
steady state. In this contribution, we partially solve this 
problem by showing how to calculate the retarded and 
advanced Green's functions, and hence the many-body 
density of states (DOS) within the dynamical mean-field 
theory approach ^9]. Since these Green's functions are a 
needed input into the reconstruction theorem, this can be 
viewed as an initial step toward a complete steady-state 
formalism. The formal development here is more gen- 
eral than the transient response formalism because 
it can be applied to any many-body lattice Hamiltonian 
that can be solved with a real time or real frequency im- 
purity solver in a dynamical mean field; here we show 
results for the Hubbard model solved with the numerical 
renormalization group (NRG) method. 

Formalism. Our focus is on the advanced and re- 
tarded Green's functions. Since the advanced Green's 
function is directly related to the retarded Green's func- 
tion via complex conjugation and an interchange of the 
two time variables, we will derive results for the retarded 
Green's function only. We use the Keldysh boundary 
condition for the nonequilibrium problem: starting our 
system in an equilibrium distribution at a constant tem- 
perature and then turning on the constant and spatially 
uniform electric field. We then let the system evolve 
forward in time until all transients have died off and 
we are left with the steady-state response. The elec- 
tric field is described by a spatially uniform vector po- 
tential in the Hamiltonian gauge, where the scalar po- 
tential vanishes E(i) = —dA(t)/cdt and we ignore all 
magnetic field effects that are present near the time the 
field is initially turned on. For a uniform field, we then 
have A{t) = —cEt since the field is turned on in the 
infinite past (but after the system has reached equilib- 
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rium at temperature . The vector potential is input 



into the Hamiltonian via the Peierls' substitution [llj, 
so that the nonequilibrium Hamiltonian is translation- 
ally invariant and can be described in momentum space. 
The momentum-dependent retarded Green's function is 
defined to be 

Gtit.t') = -z0(t-i')Tre-^^-{ck.W,cL(i')}+/2e,, 

. 

where the averages are taken with respect to the initial 
equilibrium Hamiltonian {2eq is the equilibrium parti- 
tion function), and the time evolution of the creation 
(•^kcr) ^^'^ annihilation (c^.^) operators (for electrons with 
momentum k and spin a) is in the Heisenberg picture. 
We will examine the case where the electric field lies 
along the diagonal of a hypercubic lattice in d dimen- 
sions E = E{1, 1, 1, ...) and then the limit where d oo. 

The retarded Green's function and self-energy sat- 
isfy a Dyson equation that is formally equivalent to the 
equilibrium Dyson equation except that all functions now 
depend on two time variables instead of just the time dif- 
ference 
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dt 



di'G^'At,m^{t,t')Gtit\t') 



where the noninteracting Green's function can be found 
exactly [H [H 
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Here, the bandstructure Ck satisfies [9| Ck = 
— linid-joo X^iLi <30s ki/v^, the second bandstructure 

Ek satisfies ek = — lim^j^oo i* X)iLi ^i^-^i/"^ (^^1 ener- 
gies will be measured in units of t* and we set c = 1), 
and the second line uses the Wigner coordinates of aver- 
age time T = {t + t')/2 and relative time trei = t — t'. We 
will always work with the paramagnetic solution, so we 
can drop the spin label on all functions, since the system 
will be symmetric between spin-up and spin-down, and 
there is no long-range magnetic order. 

The self-energy has no momentum dependence be- 
cause we are working in the infinite-dimensional limit; 
the perturbative result of Metzner [l4j and the Lan- 
greth rules [l^ show that the self-energy remains local in 
nonequilibrium as well. But working in the steady state 
does provide some additional simplifications to the time 
dependence which needs to be exploited to be able to 
solve the problem. To start, note that the noninteract- 
ing Green's function satisfies the gauge property which 
relates shifts in momentum to shifts in average time 



where we write the Green's function as a function of the 
Wigner coordinates. The second property is the Block 
periodicity property, which shows the system is periodic 
in the steady state 
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The gauge property implies that the local noninteract- 
ing retarded Green's function is independent of average 
time, and the Bloch periodicity property implies that the 
momentum-dependent noninteracting retarded Green's 
function is periodic in average time with the Bloch period 
2ti/E. Note that it is incorrect to assume that the steady 
state has no average time dependence. Indeed, the non- 
interacting steady-state momentum-dependent Green's 
functions explicitly depend on average time. 

The next step requires an ansatz that we cannot yet 
prove to be rigorously true. The ansatz is that the lo- 
cal retarded self-energy also has no average time depen- 
dence. The argument supporting this is that we start 
the dynamical mean-field theory algorithm with a self- 
energy equal to zero. Then the local Green's function is 
equal to the noninteracting Green's function and is av- 
erage time independent, hence the effective medium is 
average time independent and so is the impurity Green's 
function. Then the new self-energy will also be average 
time independent. Continuing in this fashion, we see that 
the DMFT algorithm will not introduce any average time 
dependence into the problem. What we cannot rule out 
is that there may exist solutions with an average time de- 
pendent self-energy; in this case we would not have any 
simple criterion to choose which one is the correct solu- 
tion to the nonequilibrium pro blem. When we analyzed 
the Falicov-Kimball model with a transient formal- 
ism, we saw that the transient solution does approach the 
steady-state solution at long times, which is an example 
showing that in at least some cases the self-energy does 
evolve to an average time independent self-energy at long 
times [l^l ■ If the self-energy has no average time depen- 
dence, then the dressed retarded Green's function satis- 
fies both the gauge property and the Bloch property, and 
in particular, the local dressed retarded Green's function 
is independent of average time. 

Using these results, we will perform a continuous 
Fourier transformation of all functions with respect to 
relative time, and a discrete Fourier series expansion 
with respect to average time (with Fourier frequen- 
cies Vn = nE, integer multiples of the Bloch fre- 
quency) [G^{T,trei) = J2nl dujGQ{vn, i^) exp{-iiy„T - 
iujtrei) /2tt\. The momentum-dependent Dyson equation 
in Eq. ([2]) becomes 



G^{Un,Uj) = Gj^"{l^n,Lu) 



(6) 
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X T.^(UJ + jVn - l'm)G^{l'n - t'm, W - il/,„), 
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which couples together the Green's functions at frequen- 
cies differing by multiples of the Bloch frequency; this 
equation has an underlying matrix structure to it that 
allows it to be solved in a straightforward fashion. In 
particular, we can restrict < tu < E when solving 
the equation, and determine the Green's function at all 

Now we describe the generalization of the iterative 
DMFT algorithm [3l for the nonequilibrium steady- 
state problem: (i) begin with a guess for the electronic 
self-energy (usually chosen to be zero); (ii) calculate the 
local retarded Green's function by solving Eq. ([6]) for each 
momentum point and summing over momentum (which 
requires a two-dimensional integration over e and e with 
a joint DOS p(e,e) = exp(— — e^)/7r [1^); (iii) extract 
the effective medium for the impurity problem from the 
Dyson equation for the impurity, using the local Green's 
function and the old self-energy; (iv) solve the impurity 
problem in the extracted effective medium to find the new 
impurity Green's function; (v) use the impurity Dyson 
equation with the new impurity Green's function and 
the old effective medium to extract the new self-energy; 
and (vi) repeat steps (ii)-(v) until the equations have 
converged. Since the impurity problem depends only on 
the relative time, we can perform a Fourier transforma- 
tion to real frequencies and use an impurity solver like 
the NRG to solve the impurity problem for the impu- 
rity self-energy. All of the nonequilibrium effects enter 
through the momentum-dependent Dyson equation and 
the self-consistency condition. Once the retarded Green's 
function is known, we compute the interacting DOS from 
p{ui) = — ImG'^^(i^„ = 0,uj)/n. Once the DOS has been 
determined, we can compute the first three moment sum 
rules and use them to test the overall accuracy of the 
computation along with the zeroth and first self-energy 
sum rule [l8| : we work at half- filling, so particle- hole sym- 
metry guarantees that the odd moments vanish. In all 
cases considered here, the zeroth moment of the DOS sat- 
isfied the sum rule to 1% or better, the second moment 
of the DOS had errors in the range from 1-25%, while 
the zeroth moment of the self-energy had errors in the 
15-30% range. All of these results are better than what 
is typically seen in equilibrium calculations. 

Results. For concreteness, we solve for the nonequi- 
librium steady-state response of the Hubbard model [19| 
in infinite dimensions; the Hubbard model involves elec- 
trons hopping between nearest-neighbor sites, with an 
on-site repulsion U . The equilibrium Hamiltonian is 

^eq = ^ ^kCko-Cko- + U C]j^.|Ci(._q-f Cpj^Cp^_q| , (7) 
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and the nonequilibrium Hamiltonian results from the 
Peierls' substitution (ck ek+Et for the steady-state 
problem). 

Details of the NRG algorithm appear in Ref. [2^. In 
most cases, we take A = 1.6 and keep 1600 states. We 



FIG. 1: Density of states for the infinite-dimensional Hub- 
bard model with E = 0.5. The panels run from systems that 
are metallic to insulating (when in equilibrium): (a) U — 0.5; 
(b) U = 2; (c) U = 4; and (d) (7 = 8. Note the change in the 
vertical axis size for the different panels. 



begin our results by showing calculations for a weak field 
case, where E — 0.5 in Fig. [1] The four panels show 
progressively larger values of the interaction strength U 
ranging from metals (in equilibrium) to Mott insulators. 
In panel (a), we have the weak coupling result with U — 
0.5. This behaves as expected, showing a broadening 
of the Wannier-Stark ladder delta functions, which are 
located at integer multiples of the Bloch frequency (0.5n 
here). As the interactions increase further, the minibands 
broaden and merge into a single band, but with a shape 
unlike that seen in equilibrium [panel (b) for U = 2], and 
then they start to form the upper and lower Hubbard 
bands [panel (c) for U = A and panel (d) for [7 = 8]; note 
that in panel (c) there is still a small peak appearing at 
the center of the density of states. The evolution with 
U can be easily understood. When the interactions are 
small enough, the driving field E primarily determines 
the DOS and we have a broadened Wannier-Stark ladder. 
As the broadening increases, the minibands merge into 
a single band. Increasing the interactions further brings 
in Mott physics which dominates the behavior and the 
DOS approaches the shape of the equilibrium DOS. 

Note that the peak at low frequency should not be con- 
fused with the quasiparticle peak that appears in equi- 
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FIG. 2: Density of states for the infinite-dimensional Hub- 
bard model with E — 4. The panels run from systems that 
are metallic to insulating (when in equilibrium): (a) U — 0.5; 
(b) (7 = 2; (c) (7 = 4; and (d) U = 8. Note the change in the 
vertical axis size for the different panels. 



librium. Indeed, the self-energy does not have a Fermi 
liquid form at all. The imaginary part does approach zero 
at w = for small U, but it does not have a quadratic 
behavior in uj, and instead looks more like a cusp. When 
the interaction strength gets large enough, the self-energy 
develops a sharp peak at w = reminiscent of the delta 
function peak for the Mott insulator in equilibrium. 

In Fig. [21 we plot the nonequilibrium DOS for the 
strong field case of E' = 4. Here the Bloch frequencies 
occur at An. In weak coupHng, the DOS corresponds to 
the Wannier-Stark ladder with the delta functions broad- 
ened and now split by U [panel (a) with U = 0.5]; most 
of the spectral weight lies around uj = 0, but there are 
still visible peaks around w = ±4. When U is increased, 
the minibands also merge as in panel (b) for [/ = 2. One 
can see the splitting of the delta functions continues to 
increase, producing structure at 4n ± 1 now. In panel 

(c) , we see a DOS for U ~ A that looks quite similar 
to the equilibrium DOS (but at a nonzero T). Panel 

(d) shows a new, and interesting effect. This case, with 
U ~ 8, has the split delta functions from the first bands 
at w = ±4 "meeting" at u; = 0. The system responds 
by creating a small weight quasiparticle-like peak in the 



DOS at Lu = 0. Although this looks like a quasiparticle 
peak, the self-energy does not behave like a Fermi liquid 
here, so the origin is different. Indeed, we find that the 
self-energy has a three-peak structure for small U — two 
broad peaks centered near uj — ±U and a narrow peak 
centered at a; = 0. As U is increased, the broad self- 
energy peaks remain at w « ±C^, but the low-frequency 
peak disappears. In this regime, the self-energy looks 
like it is developing a power-law cusp at low frequency. 
The case with U — 8 is anomalous, with the shape of the 
DOS differing significantly from what is seen at t7 = 6 
or [/ = 10. 

Conclusions. In this work, we have shown how to gen- 
eralize DMFT to nonequilibrium steady state situations, 
which can be mapped onto the same form of impurity 
problem that one uses to solve problems in equilibrium. 
We are able to explicitly determine the retarded and ad- 
vanced Green's functions, which we illustrate by exam- 
ining the Hubbard model driven by different magnitude 
electric fields. We find a rich array of behavior, including 
a broadening of the Wannier-Stark ladders, an evolution 
toward the equilibrium DOS when U is large enough, and 
a splitting of the Wannier-Stark delta functions when E 
is large. The self-energy also is anomalous, and never 
appears to illustrate behavior similar to that of a Fermi 
liquid — the nonequilibrium steady state simply behaves 
differently. 

Acknowledgments. We acknowledge useful conversa- 
tions with A. Hewson, A. -P. Jauho, V. S. Oudovenko, 
J. Serene, V. M. Turkowski, and V. Zlatic. This work 
is supported by the N. S. F. under grant number DMR- 
0705266. 



[1 
[2: 

[s: 

[4: 

[s: 
[6: 

[7: 

[8 
[9 

[10: 



[11 
[12: 



L. P. Kadanoff and G. Baym, Quantum Statistical Me- 
chanics, (W. A. Benjamin, Inc., New York, 1962). 
L. V. Keldysh, J. Exptl. Theoret. Phys. 47, 1515 (1964) 
[Sov. Phys. JETP 20, 1018 (1965)]. 

Y. Meir and N. S. Wingreen, Phys. Rev. Lett. 68, 2512 
(1992). 

P. Mehta and N. Andrei, Phys. Rev. Lett. 96, 216802 

(2006) 

F. B. Anders, eprint, arXiv:0802. 037 1\ (2008) . 

J. E. Han and R. J. Heary, Phys. Rev. Lett. 99, 236808 

(2007) 

P. Lipavsky, V. Spicka, and B. Velicky, Phys. Rev. B 34, 
6933 (1986). 

B. Velicky, A. Kalvova, and V. Spicka, J. Phys.: Confer. 
Series 35, 1 (2006). 

W. Metzner and D. VoUhardt, Phys. Rev. Lett. 62, 324 
(1989); U. Brandt and C. Mielsch, Z. Phys. B: Condens. 
Matter 75, 365 (1989). 

J. K. Freericks, V. M. Turkowski, and V. Zlatic, Phys. 

Rev. Lett. 97, 266408 (2006); J. K. Freericks, Phys. Rev. 

B 77, 075109 (2008). 

R. E. Peierls, Z. Phys. 80, 763 (1933). 

V. M. Turkowski and J. K. Freericks, Phys. Rev. B 71, 



5 



085104 (2005). 

[13] V. M. Turkowski and J. K. Freericks, in Strongly Corre- 
lated Systems: Coherence and Entanglement, edited by 
J. M. P. Carmelo, J. M. B. Lopes dos Santos, V. Rocha 
Vieira, and P. D. Sacramento (World Scientific, Singa- 
pore, 2007), pp. 187-210. 

[14] W. Metzner, Phys. Rev. B 43, 8549 (1991). 

[15] D. C. Langroth, in Linear and Nonlinear Electron Trans- 
port in Solids, edited by J. T. Devreese and V. E. van 
Doren (Plenum Press, New York and London, 1976), pp. 
3-32. 



[16] M. JarrcU, Phys. Rev. Lett. 69, 168 (1992). 

[17] V. M. Turkowski and J. K. Freericks, Phys. Rev. B 73, 

075108 (2006); Phys. Rev. B 73, 209902(E) (2006). 
[18] V. M. Turkowski and J. K. Freericks, Phys. Rev. B 77, 

096816 (2008). 

[19] J. Hubbard, Proc. R. Soc. London A276, 238 (1963). 
[20] R. Bulla, Th. Costi, Th. Pruschke, Rev. Mod. Phys. to 

appear (2008). 
[21] G. H. Wannier, Phys. Rev. 117, 432 (1960). 



